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Abstract. Anderson localization of classical waves in disordered media is a 
fundamental physical phenomenon that has attracted attention in the past three 
decades. More recently, localization of polar excitations in nanostructured metal- 
dielectric films (also known as random planar composite) has been subject of intense 
studies. Potential applications of planar composites include local near-field microscopy 
and spectroscopy. A number of previous studies have relied on the quasistatic 
approximation and a direct analogy with localization of electrons in disordered solids. 
Here I consider the localization problem without the quasistatic approximation. I 
show that localization of polar excitations is characterized by algebraic rather than by 
exponential spatial confinement. This result is also valid in two and three dimensions. 
I also show that the previously used localization criterion based on the gyration radius 
of eigenmodes is inconsistent with both exponential and algebraic localization. An 
alternative criterion based on the dipole participation number is proposed. Numerical 
demonstration of a localization-dclocalization transition is given. Finally, it is shown 
that, contrary to the previous belief, localized modes can be effectively coupled to 
running waves. 
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1. Introduction 

Anderson localization (AL) of classical waves in disordered systems is a fundamental 
physical phenomenon which takes place in the limit of strong resonant scattering when 
the "photon mean free path" becomes the order of or less than one wavelength (the 
Ioffe-Regel criterion) [1]. At a formal level, AL of electromagnetic or acoustic waves is 
similar to localization of electrons in disordered solids. There are, however, substantial 
physical differences. One such difference is that the motion of electrons can be finite. In 
contrast, classical waves can not be, in principle, indefinitely confined in a finite spatial 
region. In the case of electrons, one of the most important physical manifestations 
of AL (at zero temperature) is the conductor-insulator transition [2]. Localization 
of classical waves is manifested differently. If we consider an experiment in which 
classical waves are transmitted through a disordered slab, an analog of conductivity 
is the transmission coefficient. Assuming that the slab material is non- absorbing, the 
transmission coefficient can never turn to zero. However, if the waves are localized in 
the slab, the transmitted and reflected fields exhibit spatial variations at macroscopic 
scales (much larger than the wavelength) which are sample-specific (not self-averaging). 
Emission by localized modes in random positive-gain media is the basis for operation of 
random lasers [3]. Thus, the common feature of localized states of both electrons and 
classical waves is that the propagation can not be described by the Boltzmann transport 
equation or the diffusion approximation to the former. 

This paper is focused on AL of electromagnetic waves in a random structure which 
is distinctly different from either two-dimensional or three-dimensional random media. 
Namely, we will consider random planar composites (RPCs) [4]. The RPCs are made of 
small three-dimensional scatterers which are randomly distributed in a thin planar layer. 
Thus, the electromagnetic interaction in this system is essentially three-dimensional 
while the geometry is two-dimensional. The RPCs have attracted considerable recent 
attention due to the variety of potential applications, including surface-enhanced Raman 
spectroscopy of proteins [5-7]. The physical implication of AL of electromagnetic 
waves in RPCs can be best understood by considering an experiment in which the 
sample is excited by a near-field probe. If the electromagnetic states in the sample are 
strongly localized (at the particular electromagnetic frequency), the surface plasmon 
induced by the tip will also be localized and not spread over the entire sample. It 
must be emphasized, however, that there are other mechanisms that can lead to spatial 
exponential decay of surface plasmons. This includes decay due to absorption (Ohmic 
losses in the material and 

The localization-delocalization transition is expected to play a crucial role in the 
near-field tomographic imaging techniques of Refs. [8,9]. If the states are localized, each 
near-field measurements will be sensitive only to the local environment of the the tip, 
while in the opposite case, it would be sensitive to the structure of the sample far away 
from the tip. The relation between AL and transport of surface plasmons is discussed 
in detail in Sectional 
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The possibility and nature of AL of electromagnetic excitations in the RPCs 
have been investigated theoretically and numerically [4, 10]. While the conclusions 
given in these two references are somewhat conflicting, the respective theoretical 
approaches share some common features. Most importantly, localization of the SP 
eigenmodes was studied in the quasistatic approximation. However, AL is, essentially, 
an interference phenomenon [1]. Therefore, account of retardation is essential for its 
proper understanding. Second, a definition of localization length of a mode based on 
its "radius of gyration" was adopted in Refs. [4,10]. Here I argue that this definition, 
as well as the one based on radiative quality factor (mode lifetime) [11,12] can not be 
applied to the electromagnetic localization problem in the RPCs. 

Below, I discuss a number of important points concerning AL of classical waves, 
some of which are applicable only to the RPCs and some are more general. I also 
provide a numerical demonstration of the Anderson transition in the RPCs. To this end, 
I use a simple but physically relevant model of small spherical inclusions of diameter D 
embedded in a transparent dielectric host medium and randomly distributed in a plane 
inside an L x L box. Essentially, this is the model used in Refs. [11, 12]. However, I 
work in a different physical regime and use a different definition for localization. Most 
of the numerical examples shown below were obtained in the limit Z) < C <C A < I, 
where i is the average inter-particle distance and A is the wavelength. 

2. Physical Model 

Consider an infinite transparent dielectric host with N small identical spherical 
inclusions of diameter D randomly distributed in an L x L box in the xy-plane. 
The inclusions interact with a plane linearly-polarized electromagnetic wave with the 
wave number k = nu/c = 2n/X, where n is the host index of refraction, u is the 
electromagnetic frequency. We assume that A ^> D, but the relation between A and L 
is arbitrary. Thus, we do not use the quasistatic approximation. However, we do use the 
dipole approximation which is accurate in the limit of small density of inclusions and 
D < A. Assuming that the inclusions are non-magnetic, the electric dipole moments dj 
(i — 1, . . . , N) induced in each spherical inclusion satisfy the coupled-dipole equation 



where E inc is the amplitude of the incident wave, |k inc | = k, a is the polarizability of 
inclusion, rj is the radius-vector of the i-th inclusion and G{r^ rj) is the dyadic Green's 
function for the electric field in a homogeneous infinite host medium given by 




(1) 




(2) 
(3) 
(4) 
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Here I is the unit dyadic, = Tj — r i; = jr^l, = r^/r^- and ® denotes tensor 
product. 

The system of equation can be written in operator form as 

\d) = a(\E inc ) + W\d)) . (5) 

Here the Cartesian components of all dipole moments are given by = (i,a\d), where 
i = 1, . . . , N and a = x,y, z. The above relation defines the orthonormal basis \i, a). 

The 3iV-dimensional matrix W is complex symmetric and, hence, non-Hermitian. 
Since such matrices are not very common in physics, a brief review of their spectral 
properties is adduced. Eigenvalues of complex symmetric matrices are, generally, 
complex. The eigenvectors form a complete (but not orthonormal) basis unless the 
matrix is defective. A matrix is defective if one of its eigenvectors is quasi-null, e.g, 
its dot product with itself (without complex conjugation) is zero (see below). The 
geometric multiplicity of a defective matrix is less than its algebraic multiplicity. Non- 
degenerate symmetric matrices are all non-defective. A matrix can be defective as a 
result of random degeneracy. The probability of such event is, however, vanishingly 
small. Below, we assume that W is non-defective, which was the case in all numerical 
simulations shown below. Further, let \a) and \b) be two distinct eigenvectors of W with 
components (n\a) = a n and (n\b) = b n . The usual orthogonality condition (a\b) = is 
replaced by 

= (a\b) =0 . (6) 

n 

Note that the bilinear form in the above formula is defined without complex conjugation. 
Such forms are called quasi-scalar products and are denoted by (a\b), in contrast to the 
true scalar product (a\b). The quasi-scalar product of a vector with itself, (a\a) is, 
generally, a complex number, possibly zero. A vector whose quasi-scalar product with 
itself is zero is called quasi-null. At the same time, each eigenvector (including quasi-null 
vectors) can be normalized in the usual way, so that (a\a) = 1. 

Let w n and \ip n ) be the set of 3N eigenvalues and eigenvectors of W. We assume 
here that \ip n ) are normalized so that {ip n \ip n ) = 1. However, the quasi-scalar product 
(VVil^n) is, in general, a complex number. We can use the orthogonality rule (jOJ) to 
obtain the spectral solution to (0): 

E !^<fri*u (7 ) 

n {Vn\Vn){.Z-W n ) 

where z = 1/a. We note that this spectral solution has been obtained assuming there are 
no quasi-null eigenvectors. In the opposite case, spectral solution can not be obtained. 

For non-absorbing inclusions, Imz = —2k 3 /3. This equality can be readily obtained 
by observing that the extinction cross section a e = 47rfclma and the scattering cross 
section a s = (87r/c 4 /3)|a| 2 must be equal in the absence of absorption [13] (a = 1/z). 
Analogously, in the case of finite absorption, we have Im^; < —2k 3 /3. Consequently, 
energy conservation requires [14] that lmw n > —2k 3 /3. The eigenstates with lmw n = 



Anderson Localization of Polar Eigenmodes in Random Planar Composites 



5 



—2k 3 /3 are non-radiating f while the eigenstates with Imu> n ~ N(2k 3 /3) > are 
super-radiating. The radiative quality factor of the mode is defined as Q n = l/7 n , 
7„ = lm[w n /(2k 3 /3)] + 1. For a non-radiating state, 7„ — > and Q n — > oo. The coupling 
constant f n for the n-th mode is defined as f n = (-B inc |^ n )(V' n |Si Ilc )[|E inc | 2 (V' n |^ n )]- 1 ; 
the / n 's satisfy the sum rule J2 n fn = N. 

3. The Concept of Anderson Localized for Polar Eigenmodes 

In this Section I discuss in more detail the concept of AL of polar eigenmodes and 
the relation between localization and transport properties. In particular, a rationale is 
given for studying spatial properties of the eigenmodes which depend only on the sample 
geometry but not on the material properties of the host medium or inclusions. 

It is well known that the electromagnetic problem of two-component composites, if 
solved within the quasistatics, allows for an effective separation of material properties 
of the constituents and the geometry of the composite. This idea goes back to the 
Bergman-Milton spectral theory of composites [15] and has been used in many different 
settings. For example, dipolar (more generally, multipolar) excitations in aggregated 
spheres were studied in the 1980-ies by Fuchs, Claro and co-authors [16-18] using the 
spectral approach analogous to the Bergman-Milton spectral theory of composites. 

I have extended the quasistatic spectral theory of Refs. [15-18] to the case of samples 
which are not small compared to the wavelength, e.g., when the effects of retardation 
are important in Refs. [14,19]. In this case, similarly to the quasistatics, electromagnetic 
eigenstates of a two-component mixture or composite can be defined. These eigenstates 
turn out to be independent of the material properties of the constituents, but, unlike 
in the quasistatic limit, depend explicitly on k — ou/c. Yet, at a fixed electromagnetic 
frequency u, the spatial properties of the eigenstates can be studied irrespectively of the 
material properties of the constituents. In particular, one can argue that in extended 
systems the eigenstates can be either localized on several inclusions or delocalized over 
the whole sample, completely independently of the material properties. This was the 
point of view taken in Refs. [4, 10]. 

The connection between propagation of surface plasmon excitations in the system 
and the localization properties of the eigenstates can be established by examining the 
spectral solution for the case of local excitation (e.g., by a near-field tip) at the 
site ri . If the point of observation is located at the site r 2 (e.g., another near-field tip 
operating in the collection regime), the amplitude of the measured signal is proportional 
to the following Green's function: 



| Because, in this case, and assuming the particles are non-absorbing, the terms Imz and Imw„ in the 
denominator of (7J cancel each other. Physically, this corresponds to cancellation of radiative reaction 
due to the interference effects. See Refs. [13, 14] for more details. 




(8) 
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where 



/n(r 2 ,ri) 



{r 2 \ip n )(^n\ri) 

{^n\i>n) 



(9) 



Let us assume that at a given electromagnetic frequency u there are resonance modes, 
e.g., such mode that Re [2(0;) — w n {uS)\ ~ 0. This is, obviously, not always the case, and 
the above condition depends, in particular, on the material properties of the medium. 
However, if the resonance excitation of the system is, in principle, possible, summation 
in (jSJ) can be restricted to resonant modes. In this case, propagation of surface plasmons 
in the system is governed by the spatial dependence of the functions / n (r 2 , rx), where n 
indexes only resonant modes. 

In a strictly periodic infinite system, all eigenmodes are delocalized plane waves, 
which follows from general symmetry considerations. These delocalized modes form 
a truly continuous spectrum and can be labeled by a continuous index q, where q is 
the wave vector in the first Brillouin zone of the lattice (here the dimension of q is 
not specified). Note that the eigenmodes that belong to the continuous spectrum are 
non-normalizable in the usual sense. This means that the scalar product (ip q \ip q ) does 
not exists. Instead, the usual delta-function normalization {ip q \ipq>) = 6(q — q') must be 
used. Obviously, functions /q(r 2 ,ri) that correspond to the delocalized eigenstates are 
non-decaying as the distance between the points ri and r 2 increases. 

At this point, two important comments must be made. First, the fact that the the 
eigenmodes \ip q ) are delocalized plane waves and the corresponding functions / (? (r 2 ,r 1 ) 
are non-decaying does not necessarily imply that there can be no exponential spatial 
decay of surface plasmon excitations in the system. This will be illustrated later in 
this Section with several examples. Second, eigenvalues in a strictly periodic and non- 
chiral system are invariant with respect to the change q — > —q and, therefore, doubly 
degenerate. This does not lead to defectiveness of the operator W but the degenerate 
modes must be appropriately orthogonalized. 

If we now introduce disorder and break the translational symmetry of an infinite 
sample, the eigenmodes \ip n ) are no longer plane waves. Still, the modes can be either 
localized or delocalized. The fundamental property of localized modes is that such 
modes are square integrable in the sense that (f/v|Vv) < 00 (in fact, localized modes can 
always be normalized so that (ipnl^n) — 1) and, therefore, belong to a discrete spectrum. 
There are two consequences of this fact. First, a localized mode can be characterized by 
a "center of mass" or a point near which it is localized, which we denote by R. Second, 
localized modes can be labeled by a countable index. It is natural to use R itself as a 
label. Let all resonant modes be localized. Then the Green's function (jHJ) can be written 
as 



where summation is extended over resonant modes localized near the point R. It is 
clear from inspection of Eq. (jlUJl and the definition © that the above Green's function 




(10) 
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is spatially decaying when |r x — r 2 | — > oo, and that the rate of this decay depends on 
the spatial decay of the eigenmodes. 

However, it can be readily seen that the spatial decay of localized eigenmodes does 
not need to be exponential. This is because the essential requirement of localization 
(VVil^n) < oo can be satisfied even if the decay is algebraic. More specifically, it 
is sufficient that ?/>R.( r ) = (^H'b) decays faster than l/|r — R| d / 2 , where d is the 
dimensionality of the sample. Accordingly, the spatial decay of the Green's function 
(110)1 can be slower than exponential. In fact, I argue in this paper that electromagnetic 
eigenstates can not be localized exponentially. This is a consequence of algebraic 
decay of the free-space Green's function (J2J) in a transparent host medium. More 
specifically, if \ip n ) is an exponentially localized eigenstate, the equation W\^ n ) = w n \il) n ) 
can not be satisfied due to the fact that (r|W|^> n ) is asymptotically algebraic and 
(r\ip n ) is asymptotically exponential. The impossibility of exponential localization of 
electromagnetic eigenmodes can be also viewed as a consequence of the fact that there 
are no bound states for light. 

The above arguments apply to infinite samples. However, any numerical 
computation is restricted to finite systems. Modes of a finite system are always discrete 
and have a finite L 2 norm. The difference between the modes which remain discrete 
in the thermodynamic limit and the modes that become delocalized is then revealed 
by examining the rate at which these modes decay. The amplitudes of localized modes 
decay (away from their "center of mass" R) faster than l/|r — R^/ 2 . Consequently, a 
localized mode centered at a point R which is sufficiently far from the sample boundaries 
is insensitive to any change in the sample overall size or shape. In contrast, delocalized 
modes always remain sensitive to the exact shape of the boundaries. For example, a 
delocalized mode would change substantially if the system overall size is doubled while 
a localized mode will be insensitive to such change, as along as it is centered sufficiently 
far from the original boundaries. 

In light of the above, it is interesting to examine the applicability of a localization 
criterion which based on the eigenmode gyration radius, £ n (defined in Section El below). 
This parameter was used in Refs. [4, 10] within the quasistatic approximation. An 
eigenmode was considered to be localized if £ n <C L, where L is the sample characteristic 
size. The above inequality was claimed to be a consequence of exponential localization. 
First, we note that the condition <C L is not equivalent to (actually, is weaker than) 
exponential localization of the eigenmode. As is shown in Section Elbelow, it is sufficient 
that the eigenstate exhibits spatial decay faster than \/r d l 2+l for the above inequality to 
hold. But more importantly, application of this criterion can lead to a grossly inaccurate 
conclusion about delocalization of an eigenstate when it is, in fact, strongly localized. 
The reason for this is numerical. Indeed, even though localized eigenstates belong to 
a countable discrete spectrum, the spacings between eigenvalues w n can be arbitrarily 
small and tend to decrease with the sample size. Thus, two localized eigenstates with 
centered at different points Ri and R 2 can have very close eigenvalues w\ ~ W2- 
Numerically, these two eigenstates are quasi- degenerate. That means that any linear 
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superposition of these two eigenstates is also an eigenstate (within numerical precision 
of the computer). Namely, if IV'Ri) an d IV^) are two quasi-degenerate eigenstates, then 
the linear combinations 



are also eigenstates (within the numerical precision). But the gyration radii of \ifj-Rj), 
\ipn 2 ), on one hand, and of \ipr ± ), on the other, might be very different. Thus, if the 
first two states are localized with gyration radii £1,^2 *C L, the gyration radii of the 
eigenmodes \ip r ±) can be of the order of |Ri — R2I ~ L. I argue below that a more 
reliable numerical criterion of localization can be based on the so-called participation 
number. 

Thus, we have established that geometrical properties of electromagnetic 
eigenmodes play a fundamental role in AL. These eigenmodes are completely defined by 
the sample geometry and are independent of the material properties, as long as there 
is no more than two constituents in the medium. Therefore, following Refs. [4, 10] we 
focus on the properties of the eigenmodes and do not consider a specific material- and 
frequency-dependent model for the spectral parameter z. It should be noted, however, 
that polarization or electric field can be exponentially confined in space due to several 
reasons other than AL. The physical phenomena that lead to such confinement can 
strongly depend on material properties of the medium and should not be confused with 
AL. In the remainder of this Section, we briefly review several such phenomena. 

The first relevant example is exponential decay of evanescent waves. Evanescent 
decay can lead, for example, to exponential confinement of waves in one-dimensional 
periodic layered media (photonic crystals). However, evanescent waves are exponentially 
localized only in one selected direction and are oscillatory in any direction orthogonal 
to the former. 

The second example is decay due to dissipation. It is well known that a 
superposition of perfectly delocalized states can result in an exponentially decaying 
wave. Consider a superposition of one-dimensional waves with continuous wave-numbers 
q (e.g., in a one-dimensional periodic system): 



This equation would describe, for example, propagation of a surface plasmon along a 
linear periodic chain of nanospheres of period h and integration is extended over the 
first Brillouin zone of the lattice [20]. It can be shown exactly [21] that for q > k = u/c, 
lmw(q) = —2k 3 /3. (This is known as the light-cone condition; SP is not coupled to 
running waves because momentum of the photon can not be conserved. Due to the 
same reason, a propagating SP with q > k does not radiate.) We now make the usual 
quasi-particle pole approximation, namely 




(11) 




(12) 
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where q is, by definition, the solution to Re[z — w(q )] = 0, and find that the dipole 
moment decays in space as exp(—\x\/£), where the exponential scale is 

and 

5 = -[Imz + 2k 3 /3] (14) 

is a non-negative parameter characterizing the absorption strength of one isolated 
nanosphere. In general, it can be shown that 5 = in non-absorbing particles (whose 
dielectric function e(u) is purely real at a given frequency u). The important point is 
that we have obtained exponential decay in space, even though we have superimposed 
delocalized states (running waves exp(iqx)). And notice that the localization length in 
this case is directly defined by material properties through the relaxation constant 5. 
Numerical verification of Eq. is given in Ref . [20] . 

A third type of localization happens in the off- resonant or weak interaction limit. 
Mathematically, this takes place when z is large and we can neglect the term w(q) in 
the denominator of Ijllj) (or replace z — w(q) by a g-independent constant). This will 
result in polarization localized as sinc(x//i), or in the discrete the Kronecker 

delta-symbol. Physically, this is manifestation of the fact that surface plasmons can 
not propagate when the distance between polarizable particles is too large, or the 
electromagnetic frequency is very far from the resonance. 

4. Numerical Methods 

All numerical results shown below were obtained by direct numerical diagonalization 
of the interaction operator W. A FORTRAN code has been written to model a 
random RPC and to compute elements of the matrix W as well as its eigenvectors 
and eigenvalues. RPCs were generating by randomly placing particles inside a two- 
dimensional L x L box with the only requirement that each particle does not approach 
any of the previously placed particles closer that one diameter D. An attempt which did 
not satisfy this requirement was rejected and the process was repeated until the total 
of iV particles were placed inside the box. The box was considered to be embedded in 
infinite space; thus no periodic or other boundary conditions different from the usual 
scattering conditions at infinity were applied. 

Diagonalization (computation of eigenvectors and eigenvalues of W) was 
accomplished by utilization of the LAPACK subroutine ZGEEV. Recall that W is a 
3N x 3N complex symmetric matrix. In the case of RPC, it is also a block matrix. 
It contains an iV x N block whose eigenvectors correspond to excitations polarized 
perpendicular to the RPC and an 2iV x 2N block whose eigenvectors correspond to 
in-plane excitations. Each block can be diagonalized independently. The code was 
compiled and executed on an HP rx4640 server (1.6GHz Itanium-II cpu) with the 



dRew(q) 



dq 



a=ai\ 



(13) 
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Intel's FORTRAN compiler and MKL mathematical library. Diagonalization time (for 
serial execution) for a matrix of the size M scaled approximately as 30 x (M/1000) 3 sec. 
The relatively large computational time is a consequence of W not being Hermitian. 
Diagonalization of Hermitian matrices is much more computationally efficient. It should 
be noted that the procedure always returned M linearly-independent eigenvectors. 
There were no quasi-null eigenvectors and, correspondingly, W was not defective. 

5. Properties of Localized States 

We first discuss how to determine if a certain eigenstate \ip n ) is localized in the 
Andersen sense. Various definitions of localizations that has been used for electrons 
in disordered solids are reviewed in Ref. [22]. In the literature on localization of polar 
(electromagnetic) modes in disordered composites, two approaches have been adopted. 
The first approach is based on the assumption that the localization length £ n is the 
order of the gyration radius of the mode, (r 2 ) n — (r) 2 n , where (. . .)„ denotes a weighted 
average [4,10]. More specifically, we define for the n-th mode and for each inclusion 
located at r = rj the weight 



where a = x,y, z labels the Cartesian components of three-dimensional vectors. Note 
that the argument r^ in the expression m n (rj) is discrete. The basis \i, o~) was defined 
after Eq. (j^J). Since \i, a) form an orthonormal basis, we have J2i,a \h cr )(i, '\ — 1- We 
then recall that the eigenvectors \ip n ) are normalized so that (ip n \ip n ) = 1 and find that 
the weights satisfy the following sum rule: 



However, note that, in general, Sn m n( r «) 1> with the equality holding only in the 
quasistatic limit. This is because the basis of eigenvectors \x/j n ) is not orthonormal 
(although it is complete). Then, according to Refs. [4,10], the localization length £ n for 
the n-th eigenmode is defined as 



This definition is implicitly based on the assumption of exponential localization. 
However, exponential decay (in space) of the weights m n (rj) is impossible for classical 
waves. This follows already from the fact that the unperturbed Green's function (J2J) in 
non-absorbing transparent host media decays algebraically rather than exponentially. 
However, the exponential localization is not an absolute requirement for AL. Indeed, 
the essential feature of strongly localized states is that such states are discrete and 
can be, therefore, labeled by countable indices [23]. These indices can be associated, for 
example, with localization regions which are, by definition, countable. As a consequence, 




(15) 



(J 




(16) 




(17) 
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the localized states are normalized in the usual sense, implying that / m n (r)r d ~ l dr 
converges at the upper limit, where d is the dimensionality of embedding space and we 
have approximated summation over discrete variables r« by integration over the volume 
of the sample. Note that for the RPCs d = 2, even though the interaction is three- 
dimensional. Therefore, a state is localized if the weights decay faster than l/r d . In 
contrast, delocalized states belong to the true continuum (in an infinite system). Such 
states can not be normalized in the usual sense but instead satisfy (VVtlVv) = ~~ u )i 
where fi and v are continuous variables. Consequently, the above integral is diverging for 
delocalized states. Now consider the localization length defined by (|17|) . Without loss of 
generality, we can assume that the center of mass (the second term in (j!7p ) is zero. The 
first term converges if m n (r) decays faster than \/r d+2 and diverges otherwise. Thus, 
the requirement that ^„ C L is, in fact, much stronger than is necessary for localization 
(yet, is still weaker than the requirement of exponential localization). That is, some 
modes which are truly localized in the Anderson sense will appear to be delocalized 
according to the definition (fTTjl . 

To illustrate this point, I introduce a different localization parameter. Let 



N 

E 



ml (r ?: 



We will refer to M n as the participation number of the n-th eigenmode. It is directly 
analogous to the participation number defined as the inverse second moment of the 
probability density or the inverse fourth moment of the wave function [22]. It can be 
seen that, given the constraint (jlfij) . possible values of M n lie in the interval 1 < M n < N. 
Thus, for example, if all the weights are equal, m n (rj) = 1/N, we have M n = N. If the 
mode is localized on just one inclusion i = i , so that m n (rj) = 5j )io , we have M n = 1. 
In general, a mode can be considered localized if M n -C N. 

In Fig. ^ the participation number M n is compared to £ n for all modes in an RPC 
consisting of N = 4, 000 inclusions. While there is positive correlation between M n and 
£ n (the correlation coefficient r c is indicated in the figure), it can be readily seen that 
many modes which are localized in the sense that M n <C N have the gyration radius of 
the order of L. Thus, while some correlation between M n and £ n exists, there is virtually 
no such correlation when M n <C N. 

The second approach to defining localization which has been used in the literature 
is based on the eigenmode radiative quality factor [11, 12] Q n . Again, this definition 
rests on the analogy with bound states in quantum mechanics and the assumption of 
exponential localization. However, it is easy to see that propagating modes in three- 
dimensional transparent periodic or homogeneous media are all strictly non-radiating 
(with 7„ = 0). An example of a propagating mode in a one-dimensional periodic chain 
which is strictly non-radiating, is given in Ref. [21]. On the other hand, radiating modes 
in an RPC can, in principle, be localized. This is illustrated in Fig. El Here we plot 
the inverse radiative quality factors 7 n vs the corresponding values of M n for the same 
set of parameters as in Fig. [T] First, it can be seen that, while the localized modes 
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1 10 10 2 10 3 1 10 10 2 10 3 

Figure 1. Gyration radii £„ vs the the participation numbers M n of all modes with 
polarization perpendicular to the plane (a) and parallel to the plane (b) . Parameters: 
N = 4000, L = 1000D and A = 1QQD. Numerical value of the correlation coefficient 
r c between M n and £„, calculated for all modes, is indicated in each plot. 
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Figure 2. Inverse radiative quality factors 7„ vs M n . Same parameters as in Fig.^ 
Numerical value of the correlation coefficient r c between M n and 7„ is indicated in 
each plot. 

tend to be of higher quality, the correlation is not very strong (numerical values of the 
correlation coefficient r c are indicated in each plot). Second, there are two visibly distinct 
"branches" in Fig. Efa) (E inc perpendicular to the RPC). The top, lower quality, branch 
corresponds to modes with non-vanishing dipole moments. According the criterion 
M n <C N, a significant number of such modes is localized. 

6. Localization-Delocalization Transition 

So far, we have seen that some of the modes are localized on just a few inclusions. Now 
we investigate if these modes actually form a band. To this end, we plot in Figs. El and 
H] the values of M n vs the appropriate spectral parameter of the theory, which is the 
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real part of the corresponding eigenvalue w n . To see that Rew n is, indeed, the spectral 
parameter analogous to energy, consider the following. The n-th mode is resonantly 
excited at an electromagnetic frequency u such that Re [2 (a;) — w n (cu)] = while for 
an isolated spherical inclusion the resonance condition is Re[z(o;)] = 0. Thus, the real 
parts of the eigenvalues describe shifts of resonant frequencies of collective excitations 
relative to the respective value in the non-interacting limit. This can be illustrated with 
the following simple example. Let the polarizability of a single inclusion be given by the 
Lorenz-Lorentz formula with the first non- vanishing radiative correction [24], namely, 



where is the dielectric constant of the transparent host and e is the dielectric constant 
of the inclusions. Further, let e be given by the Drude's formula 



where T is the relaxation and eo is the intra-band input to the dielectric function. For 
simplicity, let us also assume that eo = (this will not influence any conclusions in 
a significant way). Then the resonance condition for the n-th eigenmode takes the 
following form: 



Optical resonance for an isolated spherical inclusion takes place at the Frohlich frequency 
ujp = Up/ \^3eh. The corresponding resonance mode is characterized by Rew n = 0. 
Electromagnetic interaction of the inclusions results in appearance of eigenmodes which 
are characterized by Kew n 7^ 0. Corresponding spectral resonances take place at 
frequencies different from top. Using the above model for z = 1/a and e, we can 
estimate that the spectral shifts shown in Figs. EE) are limited to OSQup < uj < l.Vlujp 
which corresponds to —2 < D 3 R,ew n < 2. Note that much larger spectral shifts can 
be obtained for larger densities of inclusions. However, consideration of larger densities 
requires that calculations are carried out beyond the dipole approximation. 

Modes polarized perpendicular to the RPC are shown in Fig. El The data for parallel 
polarization are shown in Fig. 0] Analysis of Figs. clearly reveals a transition from 
delocalized to localized states. In particular, all states with sufficiently large values of 
-D 3 |Rew; n | are localized. Such states are characterized by relatively strong interaction. In 
the case of low density (N = 10 3 , £ ~ 32D), most of the localized states are binary, i.e., 
involve excitation of only two inclusions. As the density of inclusions increases, localized 
modes involving three, four and more inclusions emerge. For eigenmodes polarized 
perpendicular to the RPC plane and in the spectral region Kew n > 0, there are also 
eigenstates with M n = 1+p, p< 1. It can be argued that such modes are localized 
on just one inclusion. Yet, Reu> n for such modes is significantly shifted from the non- 
interacting limit Kew n = This result may seem to be contradictory. Indeed, if the 
eigenmode amplitude is very small on all but just one inclusion, the latter may be seen as 




(19) 




(20) 




(21) 
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Figure 3. Participation numbers M n vs dimensionless spectral parameters D 3 Rew n 
for different densities of inclusions. Polarization is perpendicular to the RPC, other 
parameters same as in Figs. Handle 
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Figure 4. Same as in Fig [3] but for polarization being parallel to the RPC 
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not interacting with its environment. The contradiction is resolved as follows. Consider 
a mode with eigenvalue w which is localized on the i-th inclusion. The corresponding 
eigenvector must then satisfy 

w (ia\1>) = ££H^Ut)(,7t|V>> . (22) 

Here a, r label the Cartesian components of vectors and the index n that labels 
eigenmodes is omitted; we thus focus attention only on the selected eigenstate. The 
eigenmode is localized on the i-th inclusion if m(rj) = 1 — p where p < 1. Therefore 
m(rj) ~ 1 and, since the weights satisfy the sum rule (|16|). we also have m(rj) ~ p/N 
for j ^ i. We now recall that the weights are quadratic in eigenvector components. 
Consequently, (jr\ip) ~ \Jp/N for j ^ i while (ic\ip) ~ 1. It then follows from (}2"2"j) that 

w ~ N{W)Jp/N = (W)^fpN , (23) 

where (W) is the appropriate average of the interaction operator in the right-hand side 
of (|22jl and the above relation is accurate only to the order of magnitude. We thus see 
that, even if p is arbitrarily small, the spectral shift may not be small in a sufficiently 
large sample (large N). This result is not specific to the RPC's but is valid in two and 
three-dimensional disordered media as well. Of course, the value of (W) will depend 
on the dimensionality of the sample and on the density of inclusions. We note that 
(W) is the same as the factor Q introduced by Berry and Percival within the mean-field 
approximation [25]. 

The phenomenon of spectrally shifted eigenstates which are localized on just one 
inclusion is explained by constructive interference and has no counterpart for electrons in 
disordered solids. This is because the analogy between the spectral parameter Rew n and 
energy is not complete. Indeed, in the case of the classic Anderson model, an electronic 
state can be localized at an anomalously deep local potential. Such state is not influenced 
in any way by values of potentials at neighboring sites since the electron is exponentially 
localized inside the potential well. But the localization phenomenon discussed here is, 
essentially, collective and depends on the particular realization of the random sample as 
a whole. Likewise, the binary states seen in Figs. I3|4I are not necessarily binary states 
of two closely situated inclusions which interact with the rest of the sample very weakly 
(although such states are also possible; see Ref. [13] for properties of isolated dimer 
states). This is evident already from the data shown in Fig. ^ Here some of the binary 
states (with M n pa 2) are characterized by large gyration radii £ n and, therefore, are not 
localized on two closely placed inclusions. 

It can also be seen that there is a spectral region (which depends on the density of 
inclusions and on polarization of the eigenmodes) where localized and delocalized modes 
co-exists simultaneously. This spectral region is especially well pronounced in the case 
of larger densities and for polarization of eigenmodes parallel to the plane of the RPC. 
Co-existence of localized and delocalized modes with very close spectral parameters has 
been previously demonstrated for random three-dimensional fractal aggregates and was 
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referred to as inhomogeneous localization [26]. Inhomogeneous localization was also 
observed in the case of RPCs [4]. However, the findings of Refs. [4,26] were based on 
the quasistatic approximation and on the definition (fTTj) of the localization length. It 
was concluded that localization is inhomogeneous in the whole spectral range. Here we 
show that inhomogeneous localization takes place only in a transitional energy band. 
The width of this transitional band as a function of system size and density of inclusions 
needs to be further investigated, especially, without the use of dipole approximation. 

7. Coupling of Localized Models to the Far Field 

We now consider the possibility of coupling of localized modes to the external field. 
It has been previously argued that, in the quasistatic limit, strongly localized modes 
can not be effectively coupled to plane waves. Therefore, such modes were referred to 
as dark [4]. However, the dark modes become coupled to external plane waves if one 
considers first non-vanishing corrections in k, i.e., goes beyond the quasistatic limit. For 
example, in Ref. [13] it was shown that the fully antisymmetrical mode of two oscillating 
dipoles (with zero total dipole moment) can be effectively coupled to an external plane 
wave in the limit kL — > 0. This coupling is explained by a small phase shift of the 
incident wave and the high quality factor of the mode. Indeed, it can be seen from (J7|) 
that under the exact resonance condition Re(z — w n ) = 0, the excited dipole moments 
become proportional to = Q n fn, where Q n = l/^ n . Thus, even in samples that 
are small compared to the wavelength, high-quality modes with zero or vanishing dipole 
moment can be effectively excited under the resonance condition. When the sample 
size is not small compared to the wavelength, even the strict resonance condition is not 
required for effective coupling. This is illustrated in Fig. EJ Here we plot the coupling 
constants f n vs M n for the same set of parameters as in Figs.^andEl Since the coupling 
constants are normalized by the condition J2 n fn — N, a mode is coupled effectively if 
f n ~ 1. The modes with / n < 1 are coupled weakly and the modes with /„ ~ N are 
coupled strongly. Correlation between M n and f n appears to be quite weak (see the 
figure for numerical values of the correlation coefficient r c ). Most importantly, it can be 
seen that a considerable fraction of localized modes is effectively coupled to the external 
wave, although only delocalized modes can be coupled strongly. 

Perhaps, the most counter-intuitive fact about the polar eigenmodes that can be 
understood only beyond the quasistatics is that the inverse radiative quality factor 7 n 
and the coupling constant /„ are not necessarily proportional to each other. In Figs. I6J7I 
we plot /„ vs 7 n for different ratios X/L and for mode polarization perpendicular to the 
RPC (Fig. EI) and parallel to the RPC (Fig. Hi)- First, consider the case of orthogonal 
polarization. In the quasistatic limit, it can be shown [14] that (for this particular 
polarization of the eigenmodes) f n = ln- This proportionality is clearly visible in the 
case X/L = 12.8 and the correlation coefficient r c between 7„ and /„, exceeds 0.999. 
However, at smaller values of X/L, there is no strict proportionality. Thus, for example, 
in the case X/L = 0.02, the correlation coefficient is rather small (r c = 0.42). It can also 
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Figure 5. Coupling constants f n vs participation numbers M n . Coupling constants 
are computed for an external plane wave ki nc = fcx and E; nc = EqZ (a) and for 
kinc = kz and Ej nc = Eq(x + y)/V2 (b). Other parameters same as in Fig. ^ 
Numerical value of the correlation coefficient r c between M n and /„ is indicated in 
each plot. 

be seen that the modes with 7 n rs 1 have coupling constants which differ by four orders 
of magnitude and can be either weakly or strongly coupled to the far field. Likewise, 
modes that are effectively coupled to the far field (f n ~ 1) can be either weakly radiating 
(j n <C 1) or strongly radiating (j n ^> 1). 

Now consider eigenmodes which are polarized in the RPC plane (Fig. |7|). There are 
two linearly-independent in-plane polarizations of the incident wave, Ei nc . While the 
radiative factors are independent of E inc , this is not so for / n 's. For any given direction 
of E inc , there is no strict proportionality between 7 n and f n even in the quasistatic limit. 
This lack of proportionality is a consequence of polarization effects. For example, there 
might be modes which are strongly coupled to incident waves polarized along the x-axis 
but not coupled to waves polarized along the y-axis. However, the polarization effects 
can be suppressed by taking the average of f n over two linearly-independent incident 
polarizations. We denote such average as (/„) and it can be shown that in the quasistatic 
limit (f n ) = 7n/2. This is indeed confirmed by the data shown in Fig. Qualitatively, 
the data in Fig. [7| are similar to those shown in Fig. |3 with even smaller correlation 
factors. 

The finding that weakly radiating modes can be effectively coupled to propagating 
waves is counter-intuitive and even may seem to contradict conservation of energy. 
Indeed, consider excitation of a mode which is effectively coupled to the far field 
but is weakly radiating by an electromagnetic wave that is "turned on" at an initial 
moment of time t = 0. Since the mode is weakly radiating, it would not contribute 
significantly to the scattered field. Thus the incident wave would seemingly pass through 
the sample without noticeable scattering or absorption. However, at a sufficiently 
large time t = T, a steady state would be reached, in which a finite electromagnetic 
energy would be transferred to plasmonic oscillations in the mode. Since the incident 
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Figure 6. Coupling constants /„ vs inverse radiative quality factors y„ for all modes 
polarized perpendicular to the RPC with TV = 4000 inclusions and different ratios 
X/L, as indicated. In all cases, L = I000Z?. The dashed blue line corresponds to the 
quasistatic result /„ = 7„. Numerical values of the correlation coefficient r c between 
7„ and f n are indicated. 



wave was not scattered or absorbed, this contradicts energy conservation. In fact, the 
contradiction is resolved by noticing that a given mode is weekly radiating only at a 
fixed electromagnetic frequency to. (Beyond the quasistatic limit, both the eigenvectors 
\ip n ) and the eigenvalues w n are functions of uj.) But the transition process described 
above necessarily involves incident waves of different frequencies, not all of which would 
pass through the sample without scattering. This would result in non-zero extinction 
of the incident power. 



8. Summary and Discussion 

Localization of polar eigenmodes in random planar composites (RPCs) has been 
studied theoretically and numerically without the quasistatic approximation. It was 
demonstrated that the localization criteria based on exponential confinement (analogy 
with electrons in solids) can not be applied to polar excitations in disordered composites 
and in RPCs in particular. This is because localization of polar eigenmodes is algebraic 
rather than exponential. Still, eigenstates with algebraically decaying tails can be 
square-integrable and discrete, and therefore, localized in the Anderson sense. 

Note that localized eigenstates whose tails decay according to power law have 
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Figure 7. Same as in Fig. but for modes polarized in plane of the RPC and the 
coupling coefficients /„ are averaged over the two orthogonal in-plane polarizations 
of the incident wave. The dashed blue line corresponds to the quasistatic result 
(/n) = 7n/2 and (. . .) denotes polarization averaging. Numerical values of the 
correlation coefficient r c between j n and (/„) are indicated. 



been also discovered for Hamiltonians which can be represented by so-called random 
banded matrices with algebraically decaying bands [27]. Elements of such matrices ay- 
decay as \i — j\~ a . Here \i — j\ can be viewed as a liner distance in a one-dimensional 
system. It was found that for a < 1 all eigenstates are delocalized while for a > 1 
all eigenstates are algebraically localized. In the critical case a — 1, the structure of 
eigenstates was found to be multifractal. These findings were in agreement with the 
results obtained earlier for random mechanical oscillators coupled by quasistatic dipole- 
dipole interaction in Refs. [28,29]. In these references, the long-range interaction was 
accounted for perturbatively, but arbitrary dimensionality of space d was considered, 
with the conclusion that all eigenstates are delocalized when a < d. We note that 
there are several important physical differences between the model considered here 
and in Refs. [27-29]. For example, in [29] two mechanical oscillators are in resonance 
irrespectively of the distance between them if their frequencies exactly coincide. In 
the model discussed here, resonances are not mechanical but electromagnetic, with 
resonance frequency strongly depending on the geometrical arrangement of dipoles f. 

| Strictly speaking, the appropriate spectral parameter for the problem of collective electromagnetic 
excitations discussed in this paper is z(uj) given by Eq. (|19fl rather than the frequency lj itself, with 
resonances taking place at frequencies satisfying one of the equations Re[z(w) — w n ] = 0. 
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Instead of solving the problem of weakly coupled oscillators of random frequencies, 
we solve the problem of electromagnetically coupled, identical polarizable particles. 
Further, a random banded matrix of size N has N 2 mathematically independent 
elements while the interaction operator W of size N (studied in this paper) depends 
on only N mathematically independent variables rj and is not banded. However, the 
mathematical reason for algebraic rather than exponential localization appears to be 
similar in both models: the algebraic spatial decay of interaction. Theory of localization 
for systems whose Hamiltonian can be represented by power-law band random matrices 
was further developed in Refs. [30-33]. It is interesting to note that electromagnetic 
interaction in the radiation zone decays as r _1 exp(ikr) and d = 2 for the RPCs. Thus, 
the rate of algebraic decay corresponds to the regime a < d. However, the interaction is 
modulated by the exponential phase factor exp(ifcr) which makes a direct comparison 
of results problematic. 

Since localized states in an RPC have algebraically decaying tails, the previously 
used localization criterion based on the gyration radius of the mode is inapplicable. 
Consequently, an alternative approach based on the participation number has been 
used in this paper. It was shown that all electromagnetic states in the RPC whose 
resonance frequencies are shifted from those of non-interacting inclusions by a value 
larger than certain threshold are localized. The band of localized states shown in 
Figs. El and |U for sufficiently large values of |Re(iu n )| can be mapped to an interval 
of electromagnetic frequencies if the material properties of the inclusions and the host 
medium are specified. It should be noted that much stronger spectral shifts will be 
observed at higher concentrations of inclusions. Consideration of the high-density limit 
will require solving the electromagnetic problem without the dipole approximation. 
When applied to large random systems, this is a very computationally demanding 
procedure, solution to which, at least at the time being, appears to be not feasible. The 
author expects that this will not influence the localization properties of the eigenmodes. 

Finally, possibility of coupling of localized modes to the far field has been studied. 
It was shown that, contrary to the previous belief, localized modes in the RPCs can be 
effectively coupled to the far field. 
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